Condensation and interaction range in harmonic boson traps: a variational approach 
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For a gas of N bosons interacting through a two-body Morse potential a variational bound of 
the free energy of a confined system is obtained. The calculation method is based on the Feynman- 
Kac functional projected on the symmetric representation. Within the harmonic approximation a 
variational estimate of the effect of the interaction range on the existence of many-particle bound 
states, and on the N-T phase diagram is obtained. 
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> ' I. INTRODUCTION 

o : 

In this paper, we continue the investigation proposed in |jj and started in |^], on the influence of realistic interactions 
on the expression for the free energy of the interacting Bose system. The method we propose is distinct from other 
treatments in two main aspects. Firstly it allows for the treatment of finite-range interatomic potentials whereas nearly 
all theoretical studies relying on a mean-field description have focused on a two-body contact potential. Secondly, the 
exact quantum statistics, both of the condensate and of the non-condensate atoms, is treated analytically at arbitrary 
temperature. 

The experimental realization of Bose-Einstcin condensation in systems of trapped, interacting bosonic atoms |^-(^] 



[ has led to renewed theoretical efforts 



141 to understand the properties of Bose gases. The need to go beyond 
the contact potential (also discussed in |15| - |17| ) appears because the ground state energy of a contact potential 
with negative scattering length (relevant for the 7 Li) is not bounded from below. The present study of finite-range 
on : potentials, based on a variational principle resulting from the Jensen-Feynman inequality, avoids this artefact. This 
variational principle (originally formulated by Feynman ^] to treat the problem of an electron in a polarizable 
medium) is extended in section II to treat many-body systems with finite-range interactions, thereby incorporating 
the quantum statistics of the particles analytically. To investigate interaction potentials different from the contact 
potential, the T-matrix formulation used in |7],|],[l6| can also be applied. In this formulation, the limit of the contact 
potential is given by the long wavelength limit. In a sense, the present method is the 'real space' complement of the 
'momentum space' T-matrix calculation: the knowledge of the pair correlation function g(r) of the model system in 
the present approach allows to effectively study the effects of the spatial dependence of the interaction potential. 

The system which we analyze in the present paper consists of a fixed number of bosonic atoms in a parabolic 
confinement. The interatomic interaction studied in detail in this paper is a Morse potential where the parameters 
are determined by the scattering length and the neutral atomic radius. Although in the most recent experiments p~9| ] 
mixtures of gases with different spin states are examined, we consider only the spin-polarized Bose gas in the present 
analysis. The results of the path integral variational method elaborated in section II, applied to the spin-polarized 
gas of bosons interacting through a Morse potential, are reported in section III, and the specific case of lithium is 
investigated using the experimentally derived interatomic potential (from spectroscopy measurements of Abraham et 
al. pp|-p2t). The discussion of the results and a comparison with other methods are presented in section IV. 



II. STUDY OF THE MORSE POTENTIAL 



The essential property of realistic interatomic interactions is that atoms repel at short distances and attract when 
they are some distance apart. In this section, we study the effect of the Morse potential t>2(r), which has these two 
main characteristics, on a collection of bosonic atoms: 

^ = (1 _ e -(r-ro)/£)2 _ Xj (1) 

where the vector r (with length r) is the difference in position vectors of the two atoms, and ro is a parameter which 
determines the range of the potential. U is a parameter which determines the strength of the potential, and L is 
a parameter which determines the 'stiffness' of the potential near its minimum. These parameters are related to 
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experimentally observable quantities. For the Morse potential the number of bound levels Ni ev and the scattering 
length a sca t are given by |23] : 
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where the square braces denote the largest positive integer smaller than the expression between brackets. The 
remaining parameters of ([!]) are determined by a least squares method using the shape of an experimentally determined 
potential [E0h22|. 



A. The Feynman-Kac variational method 

The Feynman-Kac functional is defined as an average over a Brownian motion {R(t); t ^ 0} with a variance that is 
proportional to that of standard Brownian motion by a factor yjh/m, see e.g. p4| ]. The Brownian motion provides the 
sample paths in a 3iV-dimensional configuration space (whose elements {r 1; r2,..,rjv} are denoted by r). The initial 
and final points of these paths are incorporated in the averaging symbol E r by the index and by an indicator function 
I(R(t) — r'). Using these concepts, a propagator written as 



K(r, t; r') = E r 



I(R(t) - r') exp 



satisfies the Bloch equation for distinguishable particles 

^-K(r,t;r') = ^-W 2 K(r,t;r') 



V(R(s))ds 



-V(r)K(r,t;r') 



with 



lira K(r,t;r') = S(r-r'). 



(4) 
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Using the projection techniques borrowed from |25[ and applied to confined systems in |l|j2q|, it is easy to see that 
the partition function for N identical particles at an inverse temperature [3 — l/(fcsT) is given by 



Z(N. 



(7) 



where P denotes permutations of the particle coordinates. A summation over all elements of the permutation group 
is taken. Every permutation contributes a factor £ p which is — 1 for odd permutations of fermions and 1 in all other 
cases. If the partition function Z (N,(3) and some static correlation functions of a model system, with potential 
energy Vq, can be calculated analytically, 



Z o (N,0) = ^ JdrE r ^£ P I(R(P) - P(r))expl~J^ V (R(s))ds 



(8) 



this knowledge can be used to derive an upper bound for the free energy F — —log[Z(N,/3)]//3 ) relying on the 
Jensen-Feynman inequality: 



Z(N,f3) = ± JdrEr ^e P /(i?(/3)-i 3 W)exp|-^V (i?( S ))d S |exp|-i^ 
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(V(R(s))-V (R(s)))ds 



(V(R(s))~V (R(s))}ds 



(9) 
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In this expression, the angular brackets denote the quantum statistical expectation value 



(A(R(r))) = 



*cb> / df E ? [? ei {m - p[r] } exp { - i [ v ° 



(R(s))ds A(R(r)) 



(10) 



We consider a spin polarized gas of bosons interacting through a two-body potential v 2 such as (Q) and confined by 
an anisotropic parabolic potential. The potential energy of this system is given by 



N N N 

^ = f E + *W + n^] + E E - 

i=i i=i 1=3+1 



(11) 



with m the mass of the particles, and Tj = {xj,yj, Zj} the position of the j-th boson. The partition function, the 
density and the pair correlation function for a model system with potential energy Vq given by 



N N N 

v = E? + + +E E ?( r i- r 



(12) 



were derived analytically in refs. flpq] for the isotropic case. Substituting the real potential energy for the spin 
polarized gas of interacting bosons (fLl|), and the potential energy of the trial system ( |l^ ) in the inequality (^), one 



finds 



A" 



(13) 



AT AT 



where Fo = —log[Zo(N,/3)]//3 is the free energy of the model system with partition function Zq{N,/3), wt y : 
(fl' x ) 2 + Nn/m, and R = (1/N) X),=i r i * s l ne center-of-mass coordinate. Using the pair correlation function 
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the expectation value of the two-body potential v 2 can be rewritten 
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dr v 2 (r)g(r), 



(15) 



and hence the variational free energy can be written as 

N 

2 



3 = 1 



dr v 2 (r)g(r). 



(16) 



This is the central variational formula which we will use to find the thermodynamical properties of the spin-polarized, 
parabolically confined, interacting bosons. The essential role of the pair correlation function in the evaluation of the 
expectation value of the interaction potential is clear from (|l^) . 

The required building blocks for the variational free energy (the expectation values and the pair correlation function 
for the bosonic case), obtained previously Jl],^6| for the isotropic case, have to be extended to the case of an anisotropic 
confinement potential. This anisotropic generalization is documented in the appendix. The resulting expression for 
the variational free energy is 
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(17) 



In this expression, Zq(N,P) is the partition function of the model containing N bosons at inverse temperature 
(3=1/ (ksT). In the isotropic case, Cl' x = fl' = Cl' z = fi', and the parameters W and w are the variational parameters. 
This gives a substantial simplification, since the variation with respect to f2' can be done analytically in the case of 
isotropy with f2' = f2 as the result, i.e. the variational isotropic confinement frequency equals the isotropic confinement 
frequency of the examined system. However, for the anisotropic case, expression ( |l7j ) has to be minimized with respect 
to all four parameters Q' x , £l' y , £l' z and k to find the upper bound for the free energy. 



B. Variational free energy and condensation temperature for a Morse potential 



The variational free energy associated with the Morse potential is found using fllT]). The Morse potential appears 
from the interparticle interaction as the integral of the potential times the pair correlation function: 
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Denoting = (1 — b>)(l — b l ■ 3 )/(l — b l ), b — exp{— j3hw} and a w = y/h/mw, and using the Fourier transform of 
the pair correlation function, given in the appendix, we find for the pair correlation function in the isotropic case: 
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where Zq{N,(3) is again partition function of the model system. Details on the anisotropic pair correlation function 
are given in the appendix. In this section, we assume isotropy in order not to complicate the formulas unnecessarily. 
Expression dTl) reduces to 



N(N- 1) 
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where / is the following function of a dimensionless argument: 
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The expression for the variational free energy of this system is then given by substituting ( p0[ ) in (17). In the isotropic 
case, one finds 
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where f2 equals the experimental confinement frequency and w is the remaining variational parameter and again 
a w = y/h/mw. Details for the anisotropic model can be found in the appendix. In this inequality, the expression 

obtained in for (Xw=i r |) nas been used. 
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III. RESULTS 



The variational free energy as function of w for a system with positive scattering length differs substantially from 
the case of negative scattering length. Therefore these cases are discussed separately. 

A. Morse potentials with positive scattering length 

For Morse potentials with positive scattering length, F morse (w) has only one minimum. This minimum is located 
in < w < O and the minimum value of the free energy at low temperatures is of the order of ground state energy of 
the harmonic confinement potential. The variational free energy can be used to derive the condensation temperature 
of the Bose gas. For this purpose, the specific heat is calculated from the free energy (c = —Td 2 F/dT 2 ). For Morse 
potentials with a positive scattering length, a peak appears in the specific heat as a function of temperature, indicating 
the onset of Bose-Einstein condensation. We define the condensation temperature as the temperature at which the 
specific heat reaches its maximum. In the presence of interactions, the condensation temperature T c will differ from 
the condensation temperature Tj? of the non-interacting Bose gas in the same confinement potential. This is illustrated 
in figure 1, showing the specific heat as a function of temperature for several values of the scattering length. The 
relative shift in the condensation temperature induced by the interaction is denoted by St = {T c — T®)/T®. 

Figure 2 shows the interaction induced shift St as a function of the scattering length of the Morse potential. 
Typical parameter values for the Morse potentials used in figure 2 are U = 33.36 x 10 9 Ml, ro = 1.328 x 10~ 4 a,HO, 
L = 4.794 x 10~ 5 clho with auo = {h/mVt) 1 / 2 . Adapting the range ro of this Morse potential allows to change the 
scattering length and set it to the value one wishes to study. Examples of scattering lengths appearing in experiments 
p7| , p8[ are given in table 1. A set of Morse potentials with different scattering lengths a scat was constructed, and for 
each of these Morse potentials, the interaction induced shift St in the condensation temperature was calculated. The 
results are shown in figure 2 as the full circles. It should be noted that the condensation temperature of an ideal, 
trapped Bose gas and the condensation temperature of a Bose gas interacting through a Morse potential with zero 
scattering length coincide. The full line in figure 2 shows the predicted shift in condensation temperature for a contact 
potential, as obtained from solving the Gross-Pitaevski equation ||. The results for the contact potential and the 
Morse potential approach each other for small scattering length. One could also adapt the scattering length of the 
Morse potential by adapting the parameter U. Choosing U in stead of ro as the parameter to adjust the scattering 
length had no noticeable effect on the interaction-induced shift in the condensation temperature. In the inset, the 
optimal value of the variational parameter w is shown as a function of temperature. For a repulsive Morse potential 
(flscat/ &ho > Oj full line) we find that the size of the atom cloud is expanded since w < Q. An attractive Morse 
potential (a scat < 0, discussed below, dashed line in the inset) on the other hand will contract the gas. 

B. Morse potentials with negative scatting length 

For Morse potentials with negative scattering length, F morse (w) has two minima separated by a free energy barrier. 
There is again a minimum for w of the order of O, but a second minimum is found at a much higher frequency near 
w i=s il(aHo/ r o) 2 £1 where r reflects the range of the attractive part of the interaction and ano = y/h/mil. This 
second minimum has a free energy value of the order of —NU where U reflects the depth of the attractive part of the 
interatomic potential. The average distance between the bosons is of the order of (h/mw)' 1 ^ 2 and is thus comparable 
to the range of the interatomic potential. Hence, it is plausible that the second minimum in the free energy, which 
appears for Bose gases with negative scattering length, corresponds to a many-particle bound state. We will refer to 
this state as the 'clustered' state, and to the state corresponding to the minimum with woi the order of as the 
'gaseous' state. Since the free energy in the clustered state is lower than the free energy in the gaseous state, the 
latter is metastable with respect to a transition to the clustered state. When the scattering length is not negative, the 
minimum in the free energy at a large ('cluster') value of w is not present. This property of the free energy has been 
checked numerically for a Morse potential and has been obtained analytically for any two-step square well potential 
with a non-negative scattering length. 

For the negative scattering lengths, only the specific heat associated with the gaseous state shows a peak as a 
function of temperature, and thus we find a condensation temperature only in the gaseous state, as expected. The 
interaction induced shift in the condensation temperature is opposite to the shift for potentials with positive scattering 
length, as shown in figure 2. 
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1. Phase diagram for a Bose gas with negative scattering length 



As discussed above for the Morse potentials with negative scattering length, there are in general two minima in the 
free energy. However, when the number of particles is increased at fixed temperature, we find that the minimum in 
the variational free energy associated with the gaseous state disappears above a critical number N c of particles. The 
analytic expression for this critical number at temperature zero can be found using that Z (N-l,f3)b 3l ^/Z (N,f3) -> 1 
and Pij -> for T ->■ 0, such that 

F morse (T -> 0) = hw(N - 1) + hn + 3fi (" 2 -™ 2 ) (JV _ 1} + N{ - N -^ Ue ro/L \ e ro/L f (2a w /L) - 2f(a w /L) 
2 2 iw \/2ir L 



(23) 

Then N C (T — > 0) can be found by treating the variational equation dF/dw = as an equation in w and N and finding 
the maximal N possible as a function of w. In the gaseous state a w /ro 1 and the asymptotic form of / can be 
used, yielding 

N C {T - 0) = Vggg£» a /("^ a ) e -n,/x ( I _ M, (24 ) 
v 7 5 5 / 4 L U U eWM' v 7 



with again ano — y h/mQ,. The relation of the Morse potential parameters to the scattering length a scat of the 
potential allows to rewrite ( pi]) as N C (T — * 0) = — (\/87r/5 5 ^ 4 )afl"o/«scat ~ —0.67 ano / a scat- This theoretical value 
for N c in the limiting case of T — > was also obtained by Salasnic using a Gaussian variational wave function and 
the Rayleigh-Ritz variational principle. The advantage of the present method, however, is that unlike a variational 
approach based on a trial wave function, the current technique allows to calculate N c at any temperature. The 
temperature dependence of the critical number N C {T) is shown in figure 3. Figure 3 represents a phase diagram in 
the N, T plane for the case of negative scattering length a sca t/cLHO = —6.7 x 10~ 3 , a value typical for the experiments 
on ultra cold trapped atoms (see also table 1). This specific choice for a sca t/aHO corresponds to N c = 100 at zero 
temperature. Several regions can be distinguished in this phase diagram: a region where the metastable gaseous state 
exists and is not Bose condensed, a region where the metastable gaseous state exists as Bose condensate, and finally 
a region where only the clustered state was found by the present approach. 

Keeping the temperature T fixed, we find a number of particles Nb(T) such that Bose-Einstein condensation sets 
in for Nb(T) < N. We also find the critical number of particles N C (T) above which Bose-Einstein condensation is no 
longer possible. As long as Nb(T) < N C (T), Bose-Einstein condensation can be realized at the given temperature T 
for Nb(T) < N < N C (T). As the temperature is increased, Nb(T) — > N C (T) until at the tricritical temperature T tc 
(indicated in figure 2) one has Nb(T tc ) — N c (T tc ). For temperatures above this tricritical temperature, Bose-Einstein 
condensation is no longer possible, regardless of the number of particles in the gas. 

Keeping the number of particles N fixed, we find a condensation temperature T C (N) but also a temperature T c i(N) 
such that for temperatures lower than T c i(N) only the clustered state is found by the present approach. Thus, Bose- 
Einstein condensation for a given total number of particles N (both condensate and non-condensate) can only occur 
in a temperature range T c i(N) < T < T C (N). This means that, given a fixed total number of particles, the gas can 
only be cooled to a temperature T c i(N), and not to T — 0. Since the gas has to be cooled all the way down to T = 
in order to have a condensate fraction of 100%, and since it can only be cooled to T c i(N), the maximal condensate 
fraction of a Bose gas with negative scattering length will be less than 100% for N such that T c i(N) > 0. 

In our approach the maximum total number of bosons N C (T) in the trap at a given temperature is obtained. This 
number can be used to estimate the fraction of the atoms in the condensate using 1 — {T /T C [N C (T)]} 3 where T C (N) 
denotes the condensation temperature for N bosons. The condensate fraction can then be studied as a function of the 
temperature for the scattering length used in the preceeding calculation, a scat /aHo = —6.7 x 10~ 3 . The temperature 
dependence of this condensate fraction is shown in figure 4. In the inset, the number of atoms in the condensate, 
calculated in rcf |Q within the T-matrix approach but for a different scattering length a sca t/ o-ho = —4.64 x 10~ 4 , 
is shown for comparison. Note that the decrease in the condensate fraction is much less pronounced in the case of 
Houbiers et al. 0. This may be due to the fact that in the present approach the thermodynamical properties are 
calculated as a function of the total number of atoms in the trapped gas, both condensed and non-condensed. Houbiers 
et al., on the other hand, calculated calculated the number of condensed atoms in the grand canonical ensemble, given 
a non conserved density as a reservoir. 

As a final remark, we wish to point out that the variational method presented in section II can also be applied to a 
gas of distinguishable particles (without Bose statistics). In this manner, the effects of statistics on for example the 
clustering, can be studied. For distinguishable particles, no sums over permutations must be taken, and the expression 
for the variational free energy (in isotropic confinement) simplifies to 
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where oj = a K y/coth.(H(3w/2), and f^o is the experimental confinement frequency. The variational free energy of the 
system of distinguishable particles no longer gives rise to a peak in the specific heat as a function of temperature, 
indicating that there is no phase transition to an ordered state (the condensate). As it should be, no condensation 
temperature T c > is found. But for Morse potentials with negative scattering length, the variational free energy 
will again allow for a clustered state with respect to which the gaseous state is metastable. One may wonder whether 
the critical number of particles above which the present approach only finds the clustered state depends on the Bose 
statistics: is it the same for distinguishable particles as for bosons ? The answer as given by the present calculation, 
is no: the critical number of particles N^(T) for the distinguishable particles, shown as the dotted line in figure 3, 
differs from the critical number of particles N C (T) for the bosons except at T = 0. At a given temperature, the 
critical number above which only the clustered state persists is lower for bosons than for distinguishable particles 
under the same conditions. For temperatures approaching zero, the critical number for bosons approaches that for 
distinguishable particles. 



2. Case study: The lithium condensate 



In contrast to a standard variational approach using a trial wave function, there is no simple criterion to estimate 
how close the variational upper bound for the free energy lies to the real value of the free energy of the system. In order 
to check the robustness of the variational estimate for the free energy, we focus on a case study for lithium. Instead of 
using the Morse potential, as in thcprevious subsections, we use the real interaction potential, experimentally derived 
from spectroscopy measurements |^-|22j , and compare the results for this potential with the analytical results found 
for the Morse potential. 

This implies that we have to describe firstly how we used the real interaction potential in the calculation and 
secondly how we obtained the parameters of the Morse potential for these atoms. For distances less than 0.15 nm, 
the core-region, the real interaction potential is not determined experimentally. We introduced a constant potential 
V CO re that is adjusted to the known scattering length (2^2^]. Anticipating on the results, it should be noted that 
the introduction of a constant potential does not strongly influence the results. One of the reasons may be that 
the region of the core is small compared to the region covered by a substantial value of the pair correlation. The 
remaining experimental parameters of the trap, such as the confinement frequency, were chosen in agreement with 
the experimental setup of ]3C| l. 

The Morse potential for the lithium triplet interaction is obtained as described before. The scattering length is 
fixed to its experimental value (—27.6 nm) and the number of bound levels (eleven, |$T|) is introduced. The remaining 
parameters are determined by a least squares fit to the measured triplet potential, leading to U — 33.002 x 10 9 fifi, 
r = 1.2531 x 10" 4 a HO , L = 4.5236 x 10" 5 a HO - 

The results at zero temperature are summarized in table 2. For the variational free energy in the gaseous state, no 
differences were found in the results for the Morse potential and for the experimental potential. The variational free 
energy in the gaseous state, and the value of w associated with it, are listed in columns 2 and 3 of table 2. In the next 
two columns, the free energy and the associated w are shown for the clustered state for the Morse potential. Columns 
5 and 6 show the free energy and w for the clustered state for the experimental potential. Note that at the lowest 
particle numbers, the free energy of the clustered state is increased, and even becomes larger than the free energy of 
the gaseous state. This result deserves further attention, but in our opinion it is a typical few-body problem where 
the molecular potentials have to be taken into account to higher orders than we did here. Therefore we expect that 
the model system is no longer adequate for such a small number of particles. 

The agreement between the clustered state for the Morse potential and for the experimental potential is good, and 
increases in accuracy as the number of particles increases. Also the critical number N c remains the same. The result 
for the critical number of particles for the experimental setup of (3^] , with scattering length and trapping parameters 
such that a sca tl q.ho — 4.64 x 10~ 4 , is N c = 1443 in the present approach. This number N c represents the maximum 
number of atoms in the trap that still can undergo Bose-Einstein condensation. Thus it is an upper bound to the 
average number of atoms which is reported in |30| to lie between 600 — 1300. 
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IV. DISCUSSION AND CONCLUSION 



The method proposed in this paper allows to incorporate any two-body interaction Vi in the description of the 
thermodynamics of an interacting Bose gas. From this point of view, its application is more general than the Gross- 
Pitaevskii model which takes into account a two-body contact potential v$(r) = 4:ir5(r)h 2 a sca t/m with scattering 
length a scat , and which hence neglects the range and shape of the interatomic potential. From the expression of the 
expectation value of the two-body interaction in the Jensen-Feynman inequality, it is clear that this approximation 
will be valid as long as the range of the interatomic potential (of the order of tq for the Morse potential used above) 
is much shorter than the typical length scale £ over which the pair correlation function varies. Indeed, in that case, 
the expectation value of the interatomic potential essentially depends on the value of the pair correlation function in 
the origin and on the scattering length: 

(v 2 ) = N(N ~ 1} / dr g(r)v 2 (r) — (v 2 ) rQ<<i « ^f^- g(0) J dr v 2 (r) = 2 -^^L N {N - 1) 5 (0). (26) 

However, even for short-range interactions, the typical length scale of the pair correlation function can become 
comparable to the range of the interatomic interactions. This happens for interatomic potentials which allow for a 
clustered state (a bound state involving many particles). In that case, the details of the interaction potential become 
important and models based on the Gross-Pitaevskii equation using a contact potential lead to an unphysical result: 
the free energy in these models is not bound from below. This artefact is avoided in the present treatment by taking 
into account the range of the interatomic potential. 

Bose-Einstein condensation in a gas with negative scattering length has been achieved experimentally: for a parabol- 
ically confined gas of lithium atoms a Bose-Einstein condensate of at most 600 — 1300 bosons |3^] was observed, 
comparable to the critical number for this experiment, predicted from the theory presented here. It should be noted 
that the present analysis does not describe the dynamics of a cooling Bose gas nor does it describe the dynamics of 
the formation of the clustered state. The critical number of bosons beyond which Bose-Einstein condensation is no 
longer possible may be different in the experiment due to the dynamical effects. 

In the theoretical approach presented here, as well as in the two-fluid models [p|-[l3[ for the condensate and even in 



the Monte Carlo simulations [|32 33 1, it is tacitly assumed that the system is in thermal equilibrium. It is also clear 
that the measurements are done on systems that are not in equilibrium at all; even obtaining a steady state seems 
at present out of the experimental reach Q,^5|. This means that in discussing the theoretical results we rely on 
our intuition of how the calculated equilibrium situation is reached without conflicting with the known experimental 
facts for the non-equilibrium situation. For our approach, this means that we have to argue how it is possible that a 
metastable state with a much higher energy than the ground state can be studied experimentally, while the ground 
state itself does not even show up circumstantially. For distinguishable particles it is generally accepted that the 
life-time of a metastable state is proportional to the energy-barrier between the two states and inversely proportional 
to their energy difference. There are different scenarios possible. The first one implicates gravity: as soon as a cluster 
of a few atoms is formed it leaves the trap due to the gravitational force, thereby reducing the number of atoms 
available for condensation and increasing the average energy present in the trap, because the low-energy states can 
escape now. In this scenario one is allowed to start with a number of atoms larger than the calculated N c . After 
some time of operation, clustering will reduce the number of atoms in the trap. Once the number of atoms is lower, 
the metastable state becomes available as the first excited state of the system. If the collision rate in this state does 
not differ fundamentally from the collision rate of the previous situation, the clustering will go on, and eventually 
all particles will have left the trap. If there is a change in collision rate, or if the state is almost collision free, the 
metastable state can get the necessary life-time to be available for experimentation. For 87 Rb this change in collision- 
rate was found |36| to be due to the influence of the condensate on the loss-rate of three-body recombinations. This 
characteristic was predicted by Kagan et al. |37j who suggests also to use it as the signature of the presence of the 
condensate. The physical origin of the change in collision rate is due to a change in distribution. While in the high 
temperature phase the distribution is Gaussian (at least to a good approximation) it would remain Gaussian (with 
different parameters) if the particles were distinguishable. Projecting it on the symmetric irreducible representation 
of the permutation group as required for indistinguishable particles, the distribution takes the equilibrium of the 
condensate and its excitations into account. 

In this paper, we presented a path-integral variational method which 1) exactly treats the quantum statistics at any 
temperature and 2) allows to incorporate finite range two-body interactions in the description of the thermodynamics 
of the Bose gas. After presenting this technique, we focused on a gas of parabolically trapped bosonic atoms interacting 
though a Morse potential. The effect of this interaction on the condensation temperature was studied. For Morse 
potentials with a negative scattering length, a phase diagram was derived showing the influence of the interaction on 
the region in N, T parameter space favorable to Bose-Einstein condensation. Apart from the interaction induced shift 
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in condensation temperature, the region in N, T parameter space available for Bose-Einstein condensation is reduced 
by the presence of a clustered phase. The effect of the interaction and the Bose statistics on the critical number of 
particles above which the gaseous state was no longer found in the present treatment, was calculated. We identified 
a temperature above which Bose-Einstein condensation is no longer possible regardless of the number of particles. 
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APPENDIX 



The expectation values of one-body and two-body potential energies appearing in the Jensen-Feynman inequality 
can be expressed as a function of the density 



y^wifa) = dr vi(r)n(r) with n(r) = I 
j= i J J 



/ N 
\i=l 



(27) 



and of the pair correlation function: 

N N 
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N N 



f) 



N(N — 1) 



3 = 1 1=1 



dr v 2 {r)g(r) 



with g(r) 



1 



dk 



N(N -I) J (2tt) 



— ikr 



/ N N 

EE* 

\j=l 1=1 



(28) 
(29) 



The symbol Q(expr) — 1 if the logical expression expr is true, and zero otherwise. These quantities have been discussed 
in detail in . In this appendix, the generalization of the results of [j| to the case of anisotropic confinement will be 
given. For a spin-polarized gas of bosons in an anisotropic parabolic confinement potential with frequencies Q x , Cl v , £! 2 
and a two-body potential energy (k/4) J2j i( r j ~~ Y i) 2 i it is given by 



X(r 1 ,..,r w ;/3|r 1 ,..,r w ;0) = ^C(R)^n^ x ^ H ,^(P[r,];/3|r,;0), 

P 3 = 1 

C(R) = Ko„o BA (^R;/3|ViVR;0) 



K Wx , w , Wz ( v NIL; j3\ v NIL; 0) 



(30) 
(31) 



where the first sum runs over all particle permutations P; Wi — (Of — Nk) 1 / 2 are the renormalized frequencies; 
R = (1/N) J^j r j is t ne center of mass coordinate; and K Va ., Vy ^ Vz is the single-particle path integral propagator for a 
anisotropic harmonic oscillator with frequencies v x ,v y ,v z . Hence 
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N 



n k 



,(P[r„];/3|r„;0) 



/ N 

\i=i 



1 N f 1 

— Wi^dr* C(R)-E 



N 



,Wy ,W Z (P[r„];/3|r„;0)e 



V=l 



(32) 
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P L«=i 



(34) 



where (iV, /3) indicates the partition function in the analytical model system. Each permutation can be written 
as a cyclic decomposition with Mi cycles of length 1, M2 cycles of length 2,... The sum over permutations can be 
transformed into a sum over cyclic decompositions, as in Feynman's treatment of the ideal homogeneous Bose gas 



P {Af!,Af 2 ,..} I \ I 



IMi = N 



(35) 



The constraint in the cyclic sum appears in order to guarantee that the total number of elements in all cycles of the 
cyclic decomposition equals the total number of bosons in the gas. For example, the partition sum becomes 



z (n,p) =jh[ c(r) £ [n j^m-kA e (y: 

3=1 {M U M 2 ,..} II 1 ' J \ I 



IMi = N 



Ki = K Wx 
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(r 2 ;/3|ri;0). 



(36) 
(37) 



The restriction IMi = N is prohibitive for performing the summation directly. To lift this restriction, and perform 
the summation, the generating functions are introduced: 



S(u)=J2z (N,p)u N , 
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00 II 



ikr, U N_ 
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00 INN 
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(38) 
(39) 
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JV=0 



For these generating functions, the summations and integrations can be performed analytically as for the isotropic 
case p].p6|. Using the notations bi = exp{— /3Hwi} for i — x,y,z and j3 — l/fc^T, the generating functions are given 
by 
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(43) 



with Pi t j(b) = (1 — fc'Xl — 6' "0/(1 — b ). Two schemes have been used to retrieve the expectation values from their 
generating function. The contour-integration scheme |]38|| was applied in ]3S||, and yields 
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In principle, any value can be chosen for u. Fastest numerical convergence however is achieved for the steepest descent 
solution, when u is chosen such that N = — /3~ 2 <9(lnS(u))/<9(logu). It would be incorrect to put the expectation 
value equal to the generating function itself, in which the steepest descent value u = for the parameter u in the 
free energy generating function, is substituted. 

The second method is to collect the coefficient of tt . This results in recursion relations between the expectation 
value for N particles and the expectation values for smaller number of particles: 
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(49) 



The numerical implementation of the recursion relations gives high accuracy at the cost of compu ting time. In 
practice, they become too time-consuming for more than a few thousand particles, and expressions (M4)-(Mq) have then 
to be computed. In summary, the expectation values in the path-integral formalism were calculated by introducing 
the generating functions (|38|)-(f4"o|) for these expectation values. By constructing the generating functions, one can 
analytically perform the sum over all permutations, necessary to incorporate the quantum statistics. The expectation 
value can then be extracted from the generating function using the inversion formulas (@)-(f4q) or (47)-(E9h, 
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FIGURE CAPTIONS 



Figure 1: 

The specific heat of a parabolically trapped Bose gas of 2000 atoms is shown as a function of temperature for 
several scattering lengths. In the inset the relative shift in the condensation temperature induced by the interparticle 
interaction is shown as a function of the scattering length. The condensation temperature is determined from the 
maximum of the specific heat. The full line is derived using Gross-Pitaevskii theory 0], the filled circles are obtained 
using the approach presented here. 

Figure 2: 

The interaction induced shift of the condensation temperature is shown as a function of the interaction strength 
o-scat/a-HO f° r Morse potentials (filled circles, calculated in the present approach) and contact potentials (full line, 
calculated from the Gross-Pitaevskii equation In the inset, the variational optimal value for the parameter w/Ct 
is shown as a function of temperature, for a repulsive (full line) and an attractive (dashed line) Morse potential and 
for the non-interacting, parabolically trapped Bose gas (dotted line). 

Figure 3: 

A (TV, T)-phase diagram is shown for a gas of bosons interacting through a Morse potential with negative scattering 
length a scat /aHo = —0.0067 with ano = y/h/mQ, such that N C (T = 0) = 100. The dashed line shows the 
condensation temperature T c as a function of the number of bosons. The full line shows the critical number of bosons 
N c beyond which the gaseous state no longer exists, as a function of temperature, and the dotted line shows this 
critical number for distinguishable particles under the same conditions. The temperature T tc above which BEC is not 
possible regardless of the number of bosons, is indicated: T tc — 9.43 Ml/ks, N c (T tc ) — 1363. 

Figure 4: 

N c is the maximum number of atoms in the trap which can undergo Bose-Einstein condensation. Using this 
number, and the associated condensation temperature T C (N C ) the maximal fraction of atoms in the condensate is 
shown as a function of temperature for a gas with scattering length a scat / ano — —6.7 x 10~ 3 . In the inset, the 
maximum occupation of the condensate calculated in j7| with a T-matrix approach and for a scattering length 
o* scat i 'o-ho — —4.6 x 10 -4 is shown. 
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TABLES 



Table 1: 

Typical values for system parameters in experiments on ultra cold Bose gases. For a set of chosen alkali atoms we list 
the scattering lengths a sca t (taken from ^Tj), the frequencies of the parabolic confinement potential (taken from p8fl), 
typical number of atoms for which Bose-Einstein condensation is observed (where applicable), the effective interaction 
strength given by a sca t/ 'ano and the critical number N c (found by the present method) at zero temperature (only 
applicable for negative scattering lengths), 'n.a.' stands for 'not applicable'. 



Atom 


^scat 


n n 




a S cat/a,HO 


N c 


''Li 


-1.44 ± 0.04 nm 


1300 


144 Hz 


-4.64 x 10~ 4 


1443 


23 Na, |1,- 


-1) 4.9 ± 1.4 nm 


500000 


416 Hz 


0.0048 


n.a. 


87 Rb 


4.6 ± 1.1 nm 


4500 


187 Hz 


0.0058 


n.a. 


85 Rb 


-50 nm < a scat < 


—3 nm n.a. 


~1 Hz 


-0.0050 


134 


133 Cs 


— 13 nm 


n.a. 


34 Hz 


-0.008 7 


77 



Table 2: 

Results of the variation for the free energy at temperature zero, in units mj_,i = h = Q = 1. The first column 
shows the number of particles in the parabolic confinement. The next columns show the optimal variational value for 
the parameter w and the variational free energy for (i) the gaseous state (identical results were found for the Morse 
potential and the real potential) , (ii) the Morse potential and (iii) the experimentally derived triplet potential [pp[-p2[ . 
Comments are added in the last column. 





Gaseous 




Clustered (Morse) 


Clustered (Exact) 


comment 


N 




Fx/N 


^2, morse 


F2.morse 1 


W 2 ,lith 


F 2 ,uth/N 




2 


1.00037 


1.49981 


no minimum 




no minimum 




only 


5 


1.00091 


1.49926 


no minimum 




no minimum 




gaseous 


10 


1.00184 


1.49833 


no minimum 




no minimum 




state 


20 


1.00373 


1.49647 


no minimum 




no minimum 






21 


1.00392 


1.49628 


53389.8 


4714.57 


57555.7 


4268.57 


clustered 


22 


1.00411 


1.49610 


60304.0 


3005.24 


63920.3 


2406.39 


state is 


23 


1.00430 


1.49591 


65178.6 


1100.04 


68640.5 


357.693 


metastable 


24 


1.00449 


1.49572 


69042.2 


-948.329 


72454.2 


-1831.68 


gaseous 


100 


1.01925 


1.48140 


111411. 


-2.21442 10 s 


115733. 


-2.34133 10 5 


state is 


1000 


1.31950 


1.27743 


120331. 


-3.01629 10 6 


124993. 


-3.17466 10 6 


metastable 


1443 


2.16973 


1.11880 


120623. 


-4.39386 10 6 


125297. 


-4.62399 10 6 




1444 


no minimum 




120624. 


-4.39697 10 b 


125298. 


-4.62726 10 b 


only clustered 


1500 


no minimum 




120648. 


-4.57111 10 6 


125323. 


-4.81048 10 6 


state 



Table 3: 

In this table, a description of the different length scales used in the paper is given. 



symbol 


description 


Q*scat 


scattering length 




^/(mfl) 


a w 


y/h/ (mw) 


ar 


a w y/coth[hw/(2k B T)} 


ro 


distance at which the Morse potential 




reaches its minimum 
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Tempere, Brosens, Devreese, Lemmens ; Fig.1 




Tempere, Brosens, Devreese, Lemmens ; Fig. 2 




-0.002 0.000 0.002 0.004 0.006 0.008 0.010 

a sca/ a HO 



J.Tempere, F.Brosens, L.Lemmens, J.T.Devreese ; Fig. 3 
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Phase diagram of interacting bosons (variational) 
for a morse potential with scattering length 



a scat =-0.0067 a H0 



border between metastable 
gaseous state and BEC^^; 



(N c (T tc ),TJ 




border between metastable 
BEC and region where only 
the clustered state occurs 



"Border between metastable gaseous state 
and region where only clustered state occurs, 
for the case of distinguishable particles 



\N c (T=0) 500 1000 

Number of atoms 



1500 



Tempere, Brosens, Lemmens, Devreese ; Fig.4 
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Stoof et al. 
ref. [7] 
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